# -*- coding: utf-8 -*-
###################################################
#    This file is part of blockIMH.
#
#    blockIMH is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    blockIMH is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with blockIMH.  If not, see <http://www.gnu.org/licenses/>.
###################################################
#! /usr/bin/env python

from __future__ import division
from numpy import random, exp, sqrt, power, log, newaxis, zeros, \
        zeros_like, float32, int32, repeat, minimum, ones, mean, average, sum, \
        min, prod, nan, load, arange


def createCircPermutation(p):
    """
    Circular permutations.
    """
    permutations = zeros((p, p), dtype = int32)
    for i in range(p):
        for j in range(p):
            permutations[i,j] = (i + j) % p 
    return permutations

def createNoPermutation(p):
    """
    Same order for each of the
    parallel chain.
    """
    permutations = zeros((p, p), dtype = int32)
    for i in range(p):
        for j in range(p):
            permutations[i,j] = i 
    return permutations

def createRandomPermutation(p):
    """
    Random shuffling of the
    proposals for each parallel
    chain.
    """
    permutations = createNoPermutation(p)
    for j in range(p):
        a = permutations[:, j]
        random.shuffle(a)
        permutations[:, j] = a
    return permutations

def createHalfRandom(p):
    """
    First half of the permutations
    are random shuffling.
    Second half is in reverse order of the
    first half.
    """
    permutations = createRandomPermutation(p)
    halfp = int32(p / 2)
    for j in range(halfp):
        permutations[:, j] = permutations[:, j + halfp][::-1].copy()
    return permutations

def createStratRandom(p):
    """
    Each permutation begins with a
    different proposal, then it is
    a random shuffling on the p - 1
    last elements.
    """
    permutations = zeros((p, p), dtype = int32)
    for j in range(p):
        permutations[0, j] = j
        sequence = list(arange(0, p))
        sequence.remove(j)
        random.shuffle(sequence)
        permutations[1:p, j] = sequence
    return permutations

def createTempPermutation(p):
    """
    Same order for each of the
    parallel chain.
    """
    permutations = zeros((p, p), dtype = int32)
    for i in range(p):
        for j in range(p):
            permutations[i,j] = (((i + 1) * (j + 1)) % (p + 1))  - 1
    return permutations


## Test
#p = 4
#print "no permutation"
#print createNoPermutation(p)
#print "circular permutations"
#print createCircPermutation(p)
#print "random permutations"
#print createRandomPermutation(p)
#print "half random permutations"
#print createHalfRandom(p)
#print "stratified random permutations"
#print createStratRandom(p)
#print "francois' idea for permutations"
#print createTempPermutation(p)
#
#
#
